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ABSTRACT 

The dynamics of gamma-ray burst (GRB) jets during the afterglow phase is most reliably and accurately 
modeled using hydrodynamic simulations. All published simulations so far, however, have considered only a 
uniform external medium, while a stratified external medium is expected around long duration GRB progeni- 
tors. Here we present simulations of the dynamics of GRB jets and the resulting afterglow emission for both 
uniform and stratified external media with pext ^ for k = 0, 1,2. The simulations are performed in two di- 
mensions using the special relativistic version of the Mezcal code. Common to all calculations is the initiation 
of the GRB jet as a conical wedge of half-opening angle 0o - 0.2 whose radial profile is taken from the self- 
similar Blandford-McKee solution. The dynamics for stratified external media {k - 1, 2) are broadly similar 
to those derived for expansion into a uniform external medium {k - 0). The jet half-opening angle is observed 
to start increasing logarithmically with time (or radius) once the Lorentz factor F drops below 0^'. For larger 
k values, however, the lateral expansion is faster at early times (when F > flg') and slower at late times with 
the jet expansion becoming Newtonian and slowly approaching spherical symmetry over progressively longer 
timescales. We find that contrary to analytic expectations, there is a reasonably sharp jet break in the lightcurve 
for k = 2{sL wind-like external medium) although the shape of the break is affected more by the viewing angle 
(for 9obs < do) than by the slope of the external density profile (for < ^ < 2). Steeper density profiles (i.e. 
increasing k values) are found to produce more gradual jet breaks while larger viewing angles cause smoother 
and later appearing jet breaks. The counter-jet becomes visible as it becomes sub-relativistic, and for k = 
this results in a clear bump-like feature in the light curve. However, for larger k values the jet decelerates more 
gradually, causing only a mild flattening in the radio light curve that might be hard to discern when k -2. Late 
time radio calorimetry, which makes use of a spherical flow approximation near the non-relativistic transition, 
is likely to consistently over-estimate the true energy by up to a factor of a few for k -2, but either over-predict 
or under-predict it by a smaller factor for k 

Subject headings: gamma rays: bursts - hydrodynamics - methods: numerical - relativity 



1. INTRODUCTION 

The dynamics of gamma-ray burst (GRB) outflows de- 
pends on the density distribution of the ambient medium 
as well as on the structure of the relativistic expanding 
[ ejecta (e.g.. iMeszaros et"an 119981) . Up to the deceleration 
epoch, where most of the energy is transferred to the 
shocked external medium, the dynamics is regulated by 
I the local radial structure of the ejecta, while at later times 
(as the blastwave decelerates) it mainly depends on its 
global angular structure. In the absence of characteristic 
scales, self-similar, spheri cally symmetric solutions exist 
! jBlandford & McKee|[l976l hereafter Blandford-McKee) and 
they are widely used to interpret observational data on GRB 
afterglows. However, even the simplest departure from this 
ideal model could drastically modify the afterglow behavior. 
Anisotropies in the GRB outflow, for example, affect the 
afterglow light curve when the mean jet energy per solid 
angle within the visible region evolves significantly. As the 
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jet decelerates, the relativistic beaming weakens and the 
visible region increases. If the outflow is collimated into a 
narrow jet with reasonably sharp edges, this occurs at the 
time when the bulk Lorentz factor F equals the inverse of 
jet half-opening angle Oq. A simple analytic calculation 
using the usual scaling laws leads then to a steep ening of 
the af te rglow flux decay rate, k nown as a jet break jRhoads 
T997[ ISari. Piran & HalpernI [19991: iKumar & Panaitescu 
20001) . It is however clear from numerical studies that such 



simple scalings do not provi d e an ac curate description of 
the af terglow (iGranot et all 12001): IZhang & MacFadvenl 
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van Eerten et al.ll201 lb . Such numerical studies have so far 



been limited to the case of a uniform external density while 
the interaction of relativistic GRB jets with a non-uniform 
medium remains poorly understood. 

Motivated by this, here we study the dynamics of two- 
dimensional (2D) axially symmetric impulsive jets propagat- 
ing in a spherically symmetric stratified medium of rest-mass 
density p - Ar^^ and the resulting afterglo w emission. Since 
long duration GRBs ( iGehrels et al.l 120091) have massive star 
progenitors whose winds are expected to modify their imme- 
diate s urroundings (IChevaUer et al.H2 004^jRa mirez-Ruiz et alJ 
I2005t Ivan Marie et alJl2008t iMimic a & Gia nniosll2011h . we 
consider both steady and time varying stellar winds as possi- 
ble surrounding or external media for the GRB jet evolution. 
The case k-2 corresponds to a stellar wind for a massive star 



progenitor jChevalier & LilllOOOHPanaitescu & Kumaiil2000l: 
iRamirez-Ruiz et alJl2001l: " lWu et al.ll2005l) ' with a constant ra- 
tio of its pre-explosion mass loss rate and wind velocity 
v„, in which case p = Ar^^, where A - M„7(47rv„.). However, 
since the dependence of M„ and v,t on the time f„ before the 
stellar explosion that triggers the GRB is highly uncertain, it is 
worth considering other values of k. For example, if M,v oc f^|. 
and v„. oc then the location of a wind element at the time of 
the explosion is r = fvvV,v(f,v) oc fj+'' so that t„ oc r^^*-^^*' and 
we have M„. cc r"'^i+''\ v„, oc r''/<i+'" and p oc ^-2+(«-«/(i+A). 
For a constant wind velocity (b - 0) this gives k - 2 - a, 
which corresponds to ^ = 2 for a = (constant wind mass 
flux) and ^ = 1 for a - I (linearly increasing mass flux with 
time). 

A brief description of our numerical methods and initial 
conditions for both jet and external medium models is giv- 
ing in §|2] Detailed hydrodynamic simulations of GRB jets 
interacting with k = I, 2 stratified media are presented in §[3] 
and §111 where §[3]is devoted to the jet dynamics and the re- 
sulting afterglow emission is discussed in §|4l For complete- 
ness and comparison, the interaction with a constant-density 
medium (A: - 0) is also discu ssed, although the reader is re- 
ferred to lDe CoUe etaTI (l20Tlb for a review of the current state 
of hydrodynamical modeling with k = 0. Our conclusions are 
summarized in §|5] 

2. NUMERICAL METHODS 
2.1. Code Description and Initial Conditions 

To study the dynamics of a GRB jet propagating in a 
stratified external medium, we carry out a set of two- 
dimensional simulations using the special relativistic hydro- 
dynamic (SRHD) ve rsion of the adaptive mesh refinement 
(AMR) code Mezcal (iDe Collect al.ll2011h . The Mezcal code 
integrates the SRHD equations by using a second-order (in 
space and time, except in shocks where it reduces to first order 
in space by a minmod li miter) upwind scheme based on the 
relativistic HLL method ( [Schneider et alJ[T993h . The equa- 
tion of state ( EOS), r elating enthalpy to pressure and density, 
is take n fromlRyueT al. (2006), which approximates the exact 
ISyngd (1197 Ih EOS with an error of 0.5%. This EOS prop- 
erly recovers the correct values of the adiabatic index F in the 
ultra-relativistic (F = 4/ 3) and Newtonian (F = 5/3) regimes. 
The reader is referred to iDe Colle et al.l (1201 Ih for a detailed 
description of the code and an extensive list of numerical tests. 

For the initial conditions we use a conical wedge of half- 
opening angle 9{), within which the initial radial profiles of 
pressure, density and Lorentz factor in the post-shock region 
are taken from the spherical Blandford-McKee self-similar 
solutions for a stratified medium: 
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Two-dimensional simulations with ^ = (homogeneous 
medium), k - I and k - 2 (corresponding to a steady stellar 
wind medium) are then evolved to study the lateral expansion 
and deceleration of the jet. 
To accurately study the dynamics near the jet break time, 

an initial shock Lorentz factor of Fsh.o - V2 x 20 and an 
initial half-openingjet angle Oo - 0.2 rad are selected, so that 
Fsh.o s> The isotropic equivalent energy is taken to be 
Eiso - 10^^ erg, corresponding to a total jet energy content of 
Ejgi = Eifio(l - cos Oq) ~ 2x 10^' erg. The ambient medium is 
assumed to have a density po - Aq - 1.67 x 10""^ g cm""* (for 



the case k - 0, which corresponds to po - nonipC^ with no - 
1 cm"^*), and a pressure p = r]p()C~, with 77 = lO"'". The value 
of T] has no bearing on the outcome of the simulation as long 
as the Mach number remains large, i.e. Ai ~ ?7"''"Vsh/c » 1, 
where Vsh is the shock velocity. As the simulation continues 
to evolve well into the Newtonian regime, this condition can 
be expressed as Vjh » 3 (/y/lO""^)''^ kms"'. The density 
profiles in the cases k - 1, 2 are fixed here by assuming the 
jet break radius (in the lab frame) to b e the same for all k: 
Ri(k) ^Rj{k^O). This can be rewritten (iBlandford & McKed 
[T976h as 
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where Fj = V2/0o. 

From equation (|2]i we have A^ = Ao^*(17 - 4k)/n, so that 
the density of the ambient medium is given by 



17-4;t 
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which guarantees Rj to remain unchanged for varying k. With 
this constraint, the value of the density at the jet break radius 
p(r = Rj) differs, compared to the k - case, by factors of 
13/17 and 9/17 for k = I and k - 2, respectively. In the 
simulations presented in this paper, R^ - 9.655 x 10'^ cm, 
corresponding to a jet break time of tj = Rj/c x 372 days. 

The jet is expected to begin decelerating to non-relativistic 
speeds at 
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corresponding to Jnr ~ 970, 3800 and 1 1000 days (in the lab 
frame) for A: = 0, 1 , 2 respectively. 

The simulations with k - Q,\ employs a spherical com- 
putational domain of radial and angular size {Lr,Lg) - 
(1.1 X 10'"^ cm, 7t/2) while the simulation with k — 2 uses 
{Lr,Lg) - (2.2 X 10^^ cm, n/2). The inner boundaries are 
located at (1.8, 1.2,0.3) x lO''^ cm for k = (0, 1,2), respec- 
tively. The AMR code uses a basic grid of (100,6) cells in 
the {r,6) directions, and 15 (k - 0,1) or 16 {k - 2) lev- 
els of refinement, corresponding to a maximum resolution of 
(Ar„in, A6»„i„) = (6.71 x 10^^ 1.60 x 10"^ rad). To keep 
the resolution of the relativistic thin shell A oc f'*^*^ approx- 
imately constant, the maximum nu mber of levels of refi ne- 
ment Meveis is decreased with time dDe Colle et alJ|20TTI) as 
Meveis = max[7, A^ieveLs.o - (4 - fc) log(f/fo)/ log(2)]. The simu- 
lations are halted after 150 years. We also carried out a higher 
resolution simulation (for the k - 2 case) using a basic grid 
of (1000, 16) cells in the (r, 0) directions, and 14 levels of re- 
finement. The light curves computed from this simulation are 
very similar to those those obtained from the lower resolution 
run, implying that convergence has been achieved. 

The Mezcal code is parallelised using the "Message Passing 
Interface" (MPI) library, enabling the highest resolution sim- 
ulation to be run in about two weeks on a local supercomputer 
with 160 processors, and the low resolution in about a quarter 
of that time. 

2.2. Afterglow Radiation 
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To com pute the afterglow radia tion, we use the method de- 
scribed in iDe Colle et al.l ( 1201 ll) . As the main goal of the 
current calculations is to study the effect of the jet dynam- 
ics on the afterglow lightcurves, a simple model is employed 
to calculate the emanating radiation. It assumes synchrotron 
to be the primary emitting mechanism, while ignoring self- 
absorption and inverse Compton scatte ring. Furtherrnore, a 
simple prescription for electron cooling dPe Colle et al.l20Tl1) 
is assumed , which is similar to th e one used by lOranot et^ 
(120011) and iZhang & MacFadyenl(l2009h . 

In addition to the contributions to the afterglow radiation 
computed by post-processing the results of the hydrodynam- 
ics simulations, contributions from earlier lab frame times are 
included, corresponding to the blast-wave decelerating from 

Ti = rC^r = 1) = r,h/ V2 = 200 to Ti = 20. nemx(r/Rsh) = 
1 +2(4-fe)r^j^ (1 - r/7?sh) is a self-similar vari able which quan- 
tifies the distance from the shock front (Bla ndford & McKe^ 
119761) . These are computed using the same conical wedge 
taken out of the Blandford-McKee self-similar solution that 
is used for initializing our simulations. The mapping of the 
Blandford-McKee solution is implemented by using a high 
resolution grid, starting at the position of the shock front 
(which varies with time) and sampling the Blandford-McKee 
solution at intervals of fixed AF = 0.01. The values of the 
proper density p, internal energy density eint, 4-velocity u and 
self-similar variable x replace those coming from the simula- 
tions, and are taken from the Blandford-McKee self-similar 
solution at the relevant lab frame time. In order to calcu- 
late the contributions to the observed radiation, the mapped 
jet radial structure is subsequently integrated over all angles 
{0 < 9 < 9q\ Q < <p < 2n). This procedure provides a rea- 
sonable description of the afterglow radiation at earlier times 
and it is significantly more accurate than ignoring the contri- 
butions from lab frame times preceding the start of the simu- 
lation. 

The microphysics processes responsible for field amplifica- 
tion and particle acceleration are parametrized here by assum- 
ing that the magnetic field everywhere in the shocked region 
holds a fraction eg =0.1 of the local internal energy density 
in the flow, while the non-thermal electrons just behind the 
shock hold a fraction eg = 0.1 of the internal energy, and have 
a power-law energy distribution, Nije) ye'\ with p - 2.5. 
We also assume the source to be at a redshift of z = 1, cor- 
responding to a luminosity distance of di - 2 .05 x 10^^ cm. 
The af terglow radiation code has been tested in lDe Colle et al.l 
(1201 11) . The simulation with ^ = 0, in particular, gives after- 
gl ow light curves that are nearl y identical to those computed 
bv lZhang & MacFadvenl(l2009l) . 

3. JET DYNAMICS IN A STRATIFIED MEDIUM 

Detailed hydrodynamic simulations of the evolution of a 
GRB jet in a stratified medium with k - 0,1,2 are presented 
in Figures [T] and |2] where the density and velocity contours 
of the expanding ejecta at various times are plotted. A tran- 
sient phase caused by the sharp lateral discontinuity in the 
initial conditions is observed in all cases as the shock ex- 
pands laterally and a rarefaction front moves towards the jet 
axis. This initial phase, during which shearing instabilities are 
observed to be prominent at the contact discontinuity (sepa- 
rating the original Blandford-McKee wedge material and the 
later shocked external medium), lasts for about a dynamical 
timescale and is followed by the establishment of an egg-like 
bow shock structure that persists throughout the simulations. 
The velocity quadrivector (Figure |2]i shows strong stratifica- 
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Fig. 1 . — The temporal evolution of GRB jets in a stratified medium with 
k = 0, 1,2 (top to bottom panels, respectively). The three plotted times, 
whose exact values dependent on k, have been selected so that the Blandford- 
McKee Lorentz factor in the post-shock region Ff^ = 1) is equal to 10, 5 
and 2 (left to right). Shown are logaiithmic lab frame density cuts in cm"^. 
Calculations were done in two-dimensional spherical coordinates with the 
axes corresponding to the r~ and z~ directions in units of lO'^ cm. The 
position of the shock front con'esponding to a r(x = 1) = 5 is the same for 
all k values, consistently with the normalization used in the simulations (see 
equation lo- 
tion in the direction. The expansion velocity of the jet re- 
mains mainly radial at most angles, with a non-relativistic an- 
gular component being prominent at large angles. The sub- 
structures seen in the velocity quadrivector along the z-axis 
at late times (generated by the convergence of turbulent flow) 
carry a small fraction of the energy and have a negligible ef- 
fect on the lightcurves. 

Similar resulting bow shock structures are observed for 
^ = 0, 1 and 2. However, because the rate at which mass 
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Fig. 2. — The same evolutionary sequence depicted in Figure [T]but for tlie 
absolute value of the velocity quadrivector The superposed velocity field 
arrows are represented by a gray scale color scheme linear with respect to the 
3-velocity, with dark coiTesponding to speeds ~ c and Hghter to v <K c. 

is swept-up is larger for smaller values of k, the bow shock 
lateral expansion augments with increasing k. As clearly seen 
in Figure|3] the ratio between the bow shock width and height 
as the ejecta expand changes with k. This can be understood 
as follows. Small values of k correspond to a larger increase in 
the swept-up external mass and larger decrease in the Lorenz 
factor. For the spherical case, in particular, M{< R) oc R^^''^ 
and r oc /j-(3-*)/2^ jjjg same trend should persist for the 
non-spherical case. 

The velocity quadrivector, u - Tvjc - V/S, of the expanding 
jets are shown in Figure|4]for three different angle-integrated 
quantities: mass, energy and emissivity. The mean value of 
u is larger when weighted over the energy or emissivity than 
over the shocked rest-mass until t < 10 x tj. This clearly il- 




FiG. 3. — A comparison between the bow shock structures depicted in Figure 
[T]for k = Q (black), k = I (red) and k = 2 (blue). The two times have 
been selected so that the jet has the same Lorentz factor of 10 and 5 in all 
simulations. The evolutionary scale unit of ^ct is indicated with a black 
transverse bar. The origin of the axis is located at the right bottom corner and 
the jets main direction of propagation is toward negative x in this figure. The 
simulations are normalized with respect to ct. 
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Fig. 4. — The temporal evolution (in the lab frame) of the velocity quadrivec- 
tor, u = Vv/c = FjS, in units of the jet break time. Upper panel: The difi"erent 
lines give the evolution of u within the jet when averaged over (rest-) mass, 
energy (excluding rest mass) and over the emissivity (or contribution to the 
observed flux for a distant observer along the jet axis) at 10'^ Hz for the 
evolutionary sequences shown in Figure m The non-relativistic transition 
time, ?NR(£iso). is shown in the figure as solid vertical lines for k = 0, 1,2. 
Lower panel: A comparison between the velocity quadrivector u(t) averaged 
over energy for k = b, 1,2. The non-relativistic transition time, fNR(£iso). is 
shown as a thin vertical line with the same line-style and color as the thick 
lines for the corresponding u(t). The Blandford-McKee and Sedov-Taylor 
self-similar solutions are plotted as black thin dashed lines together with the 
corresponding ~d log u/d log t slopes. 
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lustrates, in agreement with previous an alytical and numer- 
ical results limited to the case k ^ dGranot et al.l l200Tt 
IZhang & MacFadvenl2009h . that during the relativistic phase, 
most of the shocked rest mass resides in relatively slow mate- 
rial at the edges of the jet, while most of the energy is stored 
in the fastest moving material near the head of the jet. 

As illustrated in FigureH) the Blandford-McKee and Sedov- 
Taylor self-similar solutions fail to provide an adequate de- 
scription of the jet dynamics at tj < t < fNR(£^iso) with the 
disagreement becoming less pronounced before tj and after 
fNR(£iso)- Between these two limiting cases, -dlogu/dlogt 
evolves at early and late times between the two asymptotic 
slope values, as seen in the bottom panel of Figure |4] The 
evolution of -c/log u/d log t is, however, non-monotonic as it 
first increases above (3 - k)l2 and only then decreases down 
to (3 - k)l{5 - k). This behavior is mainly caused by the 
faster decrease in F compared to a spherical flow at r > tj 
due to the lateral expansion of the jet. It also relates to the 
fact that the Blandford-McKee solution depends on Ei^o while 
the corresponding Sedov-Taylor solution uses the jet's true 
energy, Zsjet and, as a result, the ratio of u{t) for these two 

limiting cases is ~ at f = fNR(£jet) and ~ 0^^^*^"** at 

t - fNR(£iso) ~ ^^'^ *'fNR(£jet)- 

Figures |5] and |6] show the resulting R±(t), R\\(t) and Oj{t) for 
k - 0,1,2 and diff'erent recipes for estimating the transverse, 
parallel and angular size scales within the jet (e.g. when aver- 
aged over mass, energy and emissivity). For all values of k the 
early lateral spreading of the jet, which starts around t ~ tj, 
is observed to initially involve only a modest fraction of the 
total energy, with the bulk of the energy reaching angles well 
above Qo at significantly later times. 

For k - Q, previous numerical simulations and analytical 
models assuming a small lateral exp ansion for t ~ Jnr (e g- 
iGranot. Ramirez-Ruiz & Loebll2005h have shown that spher- 
ical symmetry is approached on timescales much larger than 
fNR- In particular. Figures |5] shows that the growth of is 
essentially stalled at t ~ fNR while R± continues to grow as 
the flow gradually approaches spherical symmetry. For in- 
creasing k this effect is less pronounced, since R^ continues 
to increase even after fNR(£iso), albeit more slowly. This con- 
tributes to the faster growth in Oj for lower A:-values at late 
times, contrary to the opposite situation at early times (f < tj). 
This causes GRB jets expanding into steeper density profiles 
to approach spherical symmetry at progressively later times 
as argued bv lRamirez-Ruiz & MacFadvenl(l2010l) for k — 2. 

Since the rate of lateral spreadi ng of the jet i ncreases as F 
decreases (see, e.g., equation 2 of lGrano3l2007h and T(Rj) - 
0,)' is the same for all k, then the jet lateral spreading is 
expected to increase with k for R < Rj (where r{R) de- 
creases with k for a given R), while the opposite should hold 
for R > Rj (where, for a given R, r(R) incr eases with k) . 
Such a behavior is also seen in analytic models (lGranotll2007l: 
iGranot & Piran 2011). 

Figure |5] also plots the temporal evolution of rR^/R^ x 
rOj, which is observed to approach unity at f » 
tj. This shoul d be compared with the results of semi- 
analytic models (lRhoadslll997l: iSari. Piran & Halperrilll999l: 
iKumar & Panaitescull2000l e.g.). These models predict FOj ^ 
1 at f > tj, and F to decrease rapidly with lab frame time t, 
which is not observed here. In the simulations, F decreases 
rather slowly with t (as a power-law). The jet angular size 
dj (see Figure |6]l, on the other hand, is observed to increase 
only logarithmically with t for all k until the flow becomes 
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Fig. 5. — Temporal evolution (in the lab frame) of Rj^, /fy, (rSx/Ry), 
(R\\/R±) and S = d{R^IR\\)ldt in units of the jet break time. Here Sj^ and 
R|l are the transverse (cylindrical radius) and parallel (along the ; axis) scales 
of the expanding jet, respectively. The different lines give the evolution (top 
panel) of fix defined as the transverse scale of the jet that contains 50% (solid) 
or 95% (dashed) of the total energy excluding rest mass (top panel) and the 
evolution (middle, bottom panels) of {R\\) averaged over the total energy 
excluding rest mass. 

non-relativistic. 

As shown in Figure |6] the weighted mean of Gj over 
the emissivity (and to a slightly lesser extent over the en- 
ergy) remains practically constant until f/fy ~ a few, while 
the weighted mean over the shocked rest mass is signifi- 
cantly larger, in accord with earUer results ([G ranot et al. 200ll 
iPiran & Granotll200II:IZhang & MacFadveril2009) . This indi- 
cates that, as argued before, a large fraction of the swept-up 
external rest mass is concentrated at the edges of the jet, while 
most of the energy and emission lies near the head. More- 
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Fig. 6. — The evolution of the jet half-opening angle 6, as a function of time 
for different k values. The various panels give 9j derived based on the angle- 
integrated mass, energy (excluding rest mass) and emissivity (at lO" Hz). 
Also plotted is the evolution of Bj computed as the characteristic angular 
scale containing 75% or 95% of the total energy (excluding rest mass). 
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Fig. 7. — The angular distribution of the energy content in the expanding jet 
(excluding rest mass) at / = (black line), 0.5 (red), 1 (green), 2 (dark blue), 
5 (pink), 10 (light blue), 20 (yellow), 50 (orange) years. 

over, it implies that (as discussed above and seen in the tem- 
poral evolution of 6 depicted at the bottom of Figure |5]l the 
lateral expansion at early times, t < tj, is significantly faster 
for larger values of k, while the situation is reversed at late 
times. 

Figure|2]plots the temporal evolution of the energy (exclud- 
ing rest energy) per solid angle, e = dE/dQ., as a function of 
the angle from the jet symmetry axis, for ^ = 0, 1,2. At 
f > 50 yrs the energy distribution appears nearly spherical 
for all ^s. At earlier times, a clear A:-dependence trend is ob- 
served, where the energy spreads to larger solid angles faster 
for a more stratified medium, but a correlation is less evident 
when one compares e{0) for different ^-values at the same four 
velocity u rather than the same lab frame time t. 

Abundant confirmation is provided here that the dynamics 
of GRB jets are greatly modified by the radial profile of the 
surrounding circumburst density. Most analytic formalisms 
(e.g., [Rhoads 1999) derive an exponential lateral spreading 
with lab frame time or radius at f > tj, which ultimately 
erases all information about the initial jet opening angle and 
relies solely on the true energy content of the jet: £jet. No 
exponential lateral expansion is observed in our study for 
k = 1,2, consistent with previous num erical work for ex- 
pansion in a constant density medium (Granot et al."2001t 
IZhang & MacFadvea.2009; .van Eerten & MacFadvea,201 iT . 



As illustrated in Figure |6] the evolution of the jet's angular 
scale containing a constant fraction of the total energy is log- 
arithmic and is not self-similar as it retains memory of the 
initial jet opening angle. The deviation from the expected 
self- s imilar exponential lateral expansion behavior (iGruzinovl 
l2007h might be at least partly due to u rapidly decreasing with 
the polar angle from the jet symmetry axis, so that the flow 
is no longer ultra-relativistic (m » 1) as it has been previ- 
ously assumed. Even with the expectation that su ch a self- 
simila r solution would be only very slowly attained (IGruzinovl 
|2003), the maximal Lorentz factor at the head of the jet in this 
formalism is predicted to decrease exponentially with time, 
which appears to be inconsistent with our numerical results. 

The resolution of this apparent inconsistency between an- 
alytic models and numerical simulations can be attributed to 
the modest values of 6o used in the simulations, which re- 
sult in the breakdown of the analytic models, which assume 
r » 1 and dj «; 1 soon after the jet starts spreading side- 
ways r < 6q') and be fore it can reach a phase of expo- 
nentia l lateral expansion (IWvgoda et al.l20l"lT : lGranot & PiranI 
l20Tlh . In the small region in which the analytical models 
are valid: 1 «: F < there is reaso nable agreement with 
simulation results (IWygoda et al.ll20lTl) . A g eneralization of 
these analytic models to any values of F or 0j (iGranot & PiranI 
1201 ll) shows reasonable agreement with the results of simu- 
lations from the early ultra-relativistic stage to the late New- 
tonian stage. Such generalized analytic models predict that 
if the jet is initially extremely narrow then there should still 
be an early phase of exponential lateral spreading. However, 
these models make the simplifying approximation of a uni- 
form jet, while in practice u quickly drops with 0. This causes 
a breakdown of the m » 1 assumption used to derive the 
self-similar solut ion, which is onl y slowly attained even under 
ideal conditions (iGruzino ^l2007h . 

4. AFTERGLOW LIGHTCURVES 

Figure [8] shows the emerging light curves at frequencies 
ranging from the radio to gamma-rays (v - 10^, 10", 10'^*, 
10'^, 10'^, lO'^ Hz), and corresponding spectra at differ- 
ent observed times fob.s, for k = 0, 1, 2, including the ef- 
fects of electron cooling and the contribution from a mapped 
Blandford-McKee solution (with 20 < max(F) < 200). Fig- 
ure |9] shows the light curves computed for v - 10*^, 10 
and 10'^ Hz, for the two dimensional simulation and the 
Blandford-McKee conical wedge, as in Figure[8] but illustrat- 
ing the contributions to the lightcurve arising from the various 
evolutionary stages of the blast wave, quantified here by con- 
sidering the emission from lab frame times where Fsh(0/ V2i 
given by the Blandford-McKee solution, ranges between 10 
and 20, 5 and 10, 2 and 5, 1 and 2 respectively. As expected, 
lower Lorentz factors contribute to the observed flux at later 
times. A slightly more subtle effect is that at the same ob- 
served time the flux at low frequencies comes from slightly 
later lab frame times t (corresponding to a lower Blandford- 
McKee Lorentz factor rsh(f))- This is because there is a lower 
flux contribution from the sides of the jet compared to the 
center, as reflected by the fact that the afterglow image is more 
limb-brig htened at higher frequencies and less so at lower fre- 
quencies (IGranot. Piran & Sarilll999t IGranot & L oeb 2001]; 
lGranotl2008l) 7 resulting in a smaller typical angular delay time 
(tg - R(t) a! R(p-/2c) in the arrival of photons to the observer 
(which is along the jet axis in these figures). As a result, the 
flux at the same observed time fobs is dominated by larger lab 
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Fig. 8.— Light curves at v = lO', lO", 10'^ lO'', lO", lO" Hz 
(black, red, green, blue, puiple, cyan respectively; top panels) and spectra 
at Jobs = O.I, 1, 10, 100, 1000 days (black/continuous line, orange/dotted, 
blue/dashed, purple/dash-dotted, yellow/dash-dash-dotted; bottom panels) 
for the models k = 0, 1,2 (top to bottom panels), calculated including elec- 
tron cooHng and the contiibution from a mapped Blandford-McKee solution 
(with 20 < max(r) < 200). 



frame times t. 

The spectra at different observer times are shown in Figure 
[8] For all values of k, the spectra evolves from a fast cool- 
ing (with Vc < Vin and Fy oc v'^-', v"'^^, v^''^^ for v < Vc, 
Vc < V < v,„, V > Vc respectively) to a slow cooling regime 
(with v,„ < Vc and Fv oc v^^^, v""'^'^~, v"''^~ for v < Vm, 
v,„ < V < Vc, V > v„, respectively). The characteristic fre- 
quency Vm quickly drops with time with an asymptotic slope 
of -2.9, -2.6, -2 for k - 0, 1,2 respectively, (while one ex- 
pects v„, oc f-(i5-4k)/(5-k)^ which is relatively closed to our re- 
sult), while Vc increases at late times as Vc cc t, that is, with a 
slope independent on the particular stratification of the ambi- 
ent medium (for comparison, in the Sedov-Taylor regime one 

expects Vc oc f(2<:-l)/(5-i)-)^ 

As shown in ^ a jet moving in a stratified medium (with 
k - I and k - 2) decelerates to sub-relativistic speed over 
larger distances with respect to a jet moving in an homoge- 
neous medium {k - 0). The consequences of it on the light 
curve are particularly evident at radio frequencies (Figure |9]l, 
where the contribution from mildly- and sub-relativistic ma- 
terial is negligible in the k - 2 case and dominant in the ^ = 
up to f ~ lO"* days. 

Figures [8] and |9] show a pan-chromatic dip or flattening in 
the lightcurves at around half a day for k - 0, a third of a 
day for k = 1 and significantly earlier for k - 2. This feature 
is also seen in Figure [TOl which shows the temporal index 
a = -c/log Fy/d log fobs as a function of fobs, where the earli- 
est value of a is larger than expected analytically for a spheri- 
cal flow (or for a jet viewed along its axis, before the jet break 
time). Figure|9]clearly illustrates the reason for this behavior. 
It basically occurs at the point where the dominant contribu- 
tion to the observed flux switches from the Blandford-McKee 
wedge with 20 < r.sh(f) < 200 to the simulation, which cor- 
respon ds to later lab-frame t imes. As pointed out and calcu- 
lated in lDe Colle et alj|201 ll for the spherical case, the relax- 
ation of the mapping of the analytic Blandford-McKee self- 
similar solution to the numerical solution and the finite resolu- 
tion of the simulation result in a dip in the Lorentz factor that 
is gradually recovered as the shocked region becomes wider 
and thus better resolved with time. This produces a dip in the 
lightcurve, that gradually goes away as the re solution of the 
simul ation is increased (see Figures 5, 6, 7 of iDe Colle et alj 
l20Tlh . This feature is a numerical artifact of the finite resolu- 
tion of the simulation. Similar errors in the light curves were 
also present in previous simulations for the ^ = case (e.g., 
our light curve in the case k = is nearly i dentical to that 
bv IZhang & MacFadvenI 120091 as depicted in iDe CoUe et alj 
120111) . 

A smaller contribution (although not easily quantifiable) to 
the pan-chromatic dip in the lightcurve is due to the particular 
initial conditions chosen in this paper. In fact, as the jet ini- 
tially has sharp edges (a step function in the 0-direction), once 
the simulation starts there is a relaxation period occurring in 
the lateral direction on a dynamical timescale (as a rarefaction 
wave propagates from the edge of the jet towards its center). 
This lateral transient phase triggered by the sharp-edged jet is 
also imprinted in the lightcurves around the time of the dip or 
flattening, and, contrary to the limited resolution artifacts, is 
not expected to go away as the resolution is increased. This 
artifact might be less pronounced for initial conditions that 
are smoother in the lateral direction (e.g., a jet with an initial 
Gaussian angular profile). 

Apart from this early-time, artificial feature, there is the ex- 
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Fig. 9. — Afterglow lightcurves emanating at different Lorentz factors. 
The red, green, blue, purple (dashed, dotted, dashed-dotted, dashed-dotted- 
dotted) lines are the contributions to the total light curve (in black) computed 
by using the outputs of the simulations at the lab frame times where rsh(0/ V2 
(as given by the Blandford-McKee solution) ranges between 10 and 20, 5 and 
10, 2 and 5, 1 and 2 respectively. The cyan dashed-dotted lines are the con- 
tributions from a Blandford-McKee wedge with 20 < T^hit)/ VI < 200. The 
hght curves include electron cooling. 
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Fig. 10. — "Shape of the jet break", i.e. temporal decay of light curve, given 
by (1- = -d log Fy I d log <obs as a function of fobs . at three different frequencies, 
including electron cooling. 

pected pan-chromatic jet break that is present at all frequen- 
cies above v„, and is observed between a day for A: = 2 to 
several days for A: = 0. These jet break features are discussed 
in more detail below. 

4.1. Jet Breaks 

Figure [TOlplots the shape of the jet break, i.e. the temporal 
decay index of the light curve, a = -d\og Fy/dlog fobs, as a 
function of observer time, fobs, for different observed frequen- 
cies and A:- values. We shall first discuss the pan-chromatic jet 
break features at frequencies that are above the typical syn- 
chrotron frequency at the time of the jet break, v > v,„(fobs,j)- 
As shown in Figure[TOl the temporal decay of the light curve 
becomes smoother for increasing k , as derived in analytic 
models ( iKumar & Panaitescull2000l hereafter KPOO). How- 
ever, the steepening in the lightcurve occurs within a signif- 
icantly smaller observed time period than that predicted by 
analytic models. Most of the increase in a occurs over a fac- 
tor of a; 3 - 5 in time for k = (compared to a decade in time 
predicted in KPOO) and within about one decade in time for 
k - 2 (compared to four decades in time predicted in KPOO). 
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Fig. 11. — Light curves corresponding to y = 10 , 10 , 10 Hz (from 
top to bottom panels) for different viewing angles 6obs (normahzed to the 
jet initial half-opening angle 8q) and external density profiles (k = 0, 1, 2), 
with (left panels) and without (right panels) electron cooling. The lightcurves 
con'esponding to = and k = I are multiplied by 1000 and 30, respectively. 
The lightcurves include the contribution from a mapped Blandford-McKee 
solution (with 20 < F < 200) and the numerical simulation (with 1 < F < 20). 

The relatively sharper jet break (compared to analytic expec- 
tations) in a stratified medium may permit the detection of 
such a jet break. We also note that there is an "overshoot" in 
the value of the temporal decay index a just after the jet break, 
which is more promi nent for lower Ac-values (in agreement 
with previous results; lGranotj|2007h . After this overshoot a 
gradually decreases, and there is also a noticeable curvature 
in the lightcurve as the flow becomes mildly relativistic and 
eventually approaches the Newtonian regime. The effects of 
electron cooling on the shape of the jet break appear to be 
rather modest in most cases. 

At low frequencies, v < v,„(fobs.j) (see Figure [TOl upper 
panel), there is only a very modest increase in a near fobs.j- On 
the other hand, when the break frequency v,„ sweeps past the 
observed frequency v, a very sharp break is seen (i.e. increase 
in a). Both features are present for A: = 0, and we find here 
that they also persist for higher A:-values. Moreover, we also 
find that this break is sharper for smaller A:-values. This is 
because the corresponding spectral break (at v„,) is very sharp 
for our simple broken power-law spectral emissivity model, 
and is not degraded by the contribution from multiple parts of 
the jet at smaller A:-values (in addition v,„ decreases somewhat 
faster in time at fobs > fobs.j for smaller A:-values). We expect 
that a more realistic synchrotron emissivity function would 
result in a significantly smoother spectral break at v,„, which 
would in turn lead to a correspondingly smoother temporal 
break. 

Figure [TT] shows afterglow lightcurves for three different 
observed frequencies (v - 10^, lO'-', 10'^ Hz; top to bottom 
panels), external density profiles (A: = 0, 1, 2), and viewing 
angles (^obs/^o - 0, 0.5, 1), both with and without electron 
cooling (left and right panels, respectively). Figure [12] shows 
the corresponding values of the temporal decay index a for 
V = 10'^ Hz. Figures [TTI and [T2l show that the shape of the 
jet break is predominantly regulated by the change in view- 
ing angle (within the initial jet aperture, < 0obs/^^o ^ 1) 
rather than by the external density power-law index k (in the 
range < A; < 2). For 0obs = most of the steepening occurs 
within a factor of ~ 2 - 4 in time for A: = 0, 1, 2 while for 
^obs/^o ~ 0.5 - 1 it takes ~ 1 - 2 decades for A: = 0, 1,2. This 
is particularly interesting because previous analytical work 
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Fig. 12. — "Shape of the jet break", i.e. temporal decay of light curve, given 
by a = -d\og Fy/d log tabs as a function of tabs, including electron cooling, 
atv = 10" Hz > v,„. 

have argued that the effect of varying k should be significantly 
larger It can also be seen in Figure [12] that the jet induced 
steepening starts earlier and ends later for larger A:-values and 
for larger viewing angles (or 6obs/6o values). Also, the over- 
shoot in the value of a is larger for greater A:- values or 9ohs/6o 
values. The jet break time is also observed to occurs later 
for larger viewing angles at all values of k, and varies over a 
factor of ~ 3 - 5 for < OobJOo < 1. 

The change in the jet break duration with k is due to the 
slower evolution of T with t ot R ^ ct as well as fobs for larger 
A:-values (T oc /j('--3)/2 ^ t^^^/^^^^^^''^ for a spherical flow). For 
0obs = and V > Vm(fobs,j), the jet break duration roughly cor- 
responds to the time it takes the beaming cone to grow past the 
limb-brightened outer part of the image. If crudely neglecting 
lateral spreading (since most of the emission near the jet break 
time is from within the initial jet aperture, iPiran & Granod 
1200 Ih . so that the dominant effec t is the "missi ng emission" 
from outside the edges of the jet (iGrano i l2007h ). and requir- 
ing that the beaming cone (of angle 6 < l/F around the line 
of sight) grows by a factor of fk, then this would correspond 
to a factor of ~ j(^-^i<)i(i-k) -^^ observed time. However, the 
resulting image is more limb-brighten ed for smaller ^-values 
jGranot & Loebl llOOlt iGranoj l2008l) . and, as a result, one 
might estimate fk=o ~ 1.3, fk=\ ~ 1.4, fk=2 ~ 1-5, which 
would result in factors of ~ 2, ~ 3 and ~ 5 in the observed 
time, in rough agreement with our numerical results. 

As to the effect of the viewing angle for a fixed value of k, 
the addition to the duration of the jet break relative to 0obs - 
corresponds approximately to the time it takes the edge of 
the beaming cone (l/F) to grow from Qq to 6*0 + ^obs- Thus, 
for Oobs - do this corresponds to a factor of 2 decrease in 
F, or a factor of ~ 2'^^"^*'^'''"'''^ increase in the observed time 
(i.e. factors of ~ 6, ~ 8 and ~ 16 for A: = 0, 1, and 2, re- 
spectively). This is in rough agreement with our numerical 
results. According to this simple estimate, the duration of 
the jet break for 0obs = and k - 2 should be a factor of 
~ (2/2)*'^"2<:)/(3-/r) ^ 34 ^ 81 in time, or almost two decades 
in observed time, also in agreement with the results of our 
calculations. 

4.2. Radio Calorimetry 

Figure [13] shows the radio lightcurves (at v = 10^ Hz) for 
k -0, 1 , 2 from our two-dimensional numerical simulations 
of a double-sided jet, as well as for a spherical blast wave 
with the same true energy and a double-sided cone of fixed 
half-opening angle 6q calculated from a spherical blast wave 
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Fig. 13. — Light curve at v = lO' Hz for the 2d runs (k = 0, 1,2), for 
spherical Id simulations with E = iJjci and for a cone with half-opening angle 
6q computed from spherical Id simulations with E = E\sa. The contribution 
due to the counterjet is included in the lightcurves, and it is also shown (dotted 
curves) for the 2d simulations. 

with the same isotropic equivalent energy (where Oohs - 
in the two non-spherical cases). As expected, the lightcurves 
computed from a spherical blast wave with the same isotropic 
energy and from the two dimensional simulation match rea- 
sonably well at early times. For the double-sided jet, it can be 
seen that the bump in the lightcurve near the non-relativistic 
transition time, caused because the counter-jet (whose contri- 
bution is indicated by a dashed line) becomes visible, is much 
more prominent for low values of k and becomes significantly 
more modest for larger A:-values. This effect is caused by the 
more gradual deceleration of the jet for larger A:-values (as 
the same mass of external medium is swept-up over a larger 
range in radii), which causes the counter-jet to become visible 
more gradually, resulting in a wider, lower peak flux bump. In 
particular, for k - 2\l amounts to a fairly modest and rather 
slow flattening of the lightcurve, which might be hard to dis- 
cern observationally. This might, however, not help explain 
the lack of a clear flattening or rebrightening in the late radio 
afterglow of GRB 030329 (e.g., Pihlstrom et al. 2007), since 
in that case detailed a fterglow modeling favors a uniform ex- 
ternal density (k - 0: lvan der Horst et al.ll2008h . 

Comparison of the radio flux at late times from a double- 
sided jet and from a spherical blast wave with the same true 
energy near the non-relativistic transition time shows that they 
are broadly similar but may differ by up to a factor < 3. For 
k - Q and k - \ the spherical analog slightly over-predicts 
the flux before the contribution from the counterjet becomes 
important, and under-predicts the flux once the emission from 
the counter-jet becomes dominant, while for k — 2 the spher- 
ical analogue consistently under-predicts the flux, by up to a 
factor of < 3. This may result in an small but not negligi- 
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ble error in the estimation of the true energy in the double- 
side jet assuming a spherical sub-relativi stic flow, as is com- 
monly done in radio calorimetry studies ("Kaneko et al."2007|; 
Berger et alJ I2O0I [Frail et all 12005 : Gorosabel et al. 200^ 
Kulkarni et al.lll99 8l). both over- or under- estimating the real 
true energy depending on the stratification of the ambient 
medium and the observer time). 

5. DISCUSSION 

We have studied the dynamics of GRB jets during the af- 
terglow stage as they propagate into different external density 
profiles, pext = Ar for k - 0, 1,2, using detailed hydrody- 
namic simulations. Our main results, which relate both to the 
dynamics and the resulting afterglow emission, can be sum- 
marized as follows. 

For the same initial half-opening angle 60 and external den- 
sity at the jet break radius (which is defined by TiiRj) = 0^'), 
the lateral spreading is initially (at R < Rj) larger for higher 
^-values. This arises because at the same radius (or lab 
frame time) the typical Lorentz factor is lower At late times 
(R > Rj) the situation is reversed, and the eff'ective jet opening 
angle at a fixed lab frame time is similar for different A:-values. 
Since for higher fc-values a larger range of radii is required in 
order to sweep-up the same amount of mass, the whole evo- 
lution extends over a much wider range of radii and times. As 
a result, the jet break in the afterglow lightcurve is smoother 
and more gradual, the non-relativistic transition occurs later, 
and the flow approaches spherical symmetry more slowly and 
over longer timescales. The effective jet opening angle is ob- 
served to increase only logarithmically with lab frame time 
(or radius) once the jet comes into lateral causal contact (i.e. 
when r drops below 0y '). 

As long as the jet is relativistic, most of the energy and 
emission are concentrated near the head of the jet while the 
slower material at the edges carries relatively little energy 
(even though it carries a substantial fraction of the swept-up 
rest-mass). This holds true for all fc-values. Once the jet be- 
comes sub-relativistic, at f > rNR(£^iso), it quickly spreads lat- 
erally and swiftly starts to approach spherical symmetry. The 
energy weighted mean value of u(t) is observed to be of order 
unity at t/tj ~ 2 rather than at t ~ ?NR(£iso), as one might 
naively expect. We find that there is little A:-dependence on 
the temporal evolution of Oj, so that irrespective of the exter- 
nal medium radial profile, all of the expanding jets approach 
spherical symmetry at similar times (~ 1-1.5 decades after 



tj). A similar conclusion can be reached from the calculated 
evolution of R\\/R± with t/tj. 

We find that contrary to the expectations of analytic mod- 
els, the shape of the jet break is affected more by the viewing 
angle (within the initial jet aperture, < 0obs/^^o ^1) than by 
the steepness of the external density profile (for < k < 2). 
Larger viewing angles result in a later jet break time and a 
smoother jet break, extending over a wide range in time, and 
with a larger overshoot (initial increase in the temporal decay 
index a beyond its asymptotic value), which is observed to 
be more prominent for lower ^-values. Larger ^-values result 
in more gradual jet breaks, but the sharpness of the jet break 
is affected even slightly more by the viewing angle as argued 
above. The counter-jet becomes visible around Jnr, and for 
^ = this results in a clear bump in the light curve. However, 
for larger A:-values the jet deceleration is more gradual and as 
a result a wider and lower bump is produced, which becomes 
hard to detect for k - 2, where it reduces to a mild flatten- 
ing in the light curve. This may explain the lack of a clear 
counter-jet signature in some late time radio afterglow light 
curves of long duration GRBs although the dynamical com- 
plexity of their surrounding circumburst medium seriously 
limits th e validity of a non-evolvi ng power-law density pro- 
file (e.g. lRamirez-Ruiz et al.ll2005l) . 

Finally, we showed that the use of a spherical blast wave 
for estimating the total energy of the jet, as is commonly done 
in radio calorimetry studies, results in an error in the estima- 
tion of the true energy content of the jet that depends on the 
stratification of the ambient medium (being on average larger 
for k - 2). In particular, in the case k - 2, the spherical 
blast wave analogy consistently overestimates the true energy, 
while for the cases k - Q and k - \ \l produces and under- or 
an over-estimate depending on whether the estimation of the 
jet energy is done before or after the non-relativistic transition 
time. 
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